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Abstract 



We present an efficient algorithm for obtaining the gauge-invariant gra- 
dient expansion of the local density of states and the free energy of a clean 
superconductor. Our method is based on a new mapping of the semiclassical 
linearized Gorkov equations onto a pseudo-Schrodinger equation for a three- 
component wave- function ip{x), where one component is directly related to 
the local density of states. Because ip{x) satisfies a linear equation of motion, 
successive terms in the gradient expansion can be obtained by simple linear 
iteration. Our method works equally well for real and complex order param- 
eter, and in the presence of arbitrary external fields. We confirm a recent 
calculation of the fourth order correction to the free energy by Kosztin, Kos, 
Stone and Leggett [Phys. Rev. B 58, 9365 (1998)], who obtained a discrep- 
ancy with an earlier result by Tewordt [Z. Phys. 180, 385 (1964)]. We also 
give the fourth order correction to the local density of states, which has not 
been published before. 

PACS numbers: 74.20.-z, 74.25.Bt 
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I. INTRODUCTION 

The phenomenological Ginzburg-Landau theory has proven to be a powerful method 
to study the physical properties of superconductors in spatially varying external fields [|l|. 
The form of the free energy F{A(r)} as a functional of the complex superconducting order 
parameter A(r) follows from general symmetry arguments. Alternatively, for temperatures 
in the vicinity of the critical temperature and for slowly varying external fields, F{A{r)} 
can be derived microscopically from the Gorkov equations of superconductivity 0. The 
microscopic approach can also be used to obtain corrections to the Ginzburg-Landau free 
energy functional. Tewordt explicitly calculated F{A(r)} up to fourth order in the 
gradients of A(r). By comparing the terms with four and two gradients, he could estimate the 
range of validity of the usual approximation for the Ginzburg-Landau free energy functional, 
where only terms with two gradients are retained. Unfortunately, a direct expansion of the 
free energy of a superconductor in powers of gradients of A(r) is quite laborious, so that 
the result of Tewordt is rather difficult to verify, and the calculation of even higher orders 
seems almost impossible. 

Recently Kosztin, Kos, Stone and Leggett 0], and Kos and Stone (KKSL) devel- 
oped new and more efficient algorithms to obtain the gradient expansion of the free energy 
F{A(r)}. Starting point of their calculations are the Bogoliubov - de Gennes equations 
for the two-component wave-functions and eigen-energies of a superconductor in an exter- 
nal field. Assuming that |A(r)| is small compared with the Fermi energy, the low-energy 
physics can be obtained by linearizing the energy dispersion in the Bogoliubov - de Gennes 
equations. Then one arrives at the so-called semi-classical Andreev equations [0, which 
involve only first order derivatives. In the work ||^ KKSL then use the special properties of 
the Andreev equations to show that for real A(r) the gradient expansion of F{A(r)} can 
be indirectly obtained from the gradient expansion of the resolvent of a one-dimensional 
Schrodinger operator, which is nothing but the square of the Andreev-Hamiltonian. This 
resolvent satisfies a non-linear second order differential equation, the so-called Gelfand-Dikii 
equation 0,|^J^. KKSL proceed by solving the Gelfand-Dikii equation iteratively, and then 
use the approximate solution to reconstruct the gradient expansion of F{A{r)}. Although 
this strategy is more efficient than the direct approach used by Tewordt [^, this procedure 
still has the disadvantage that, in order to obtain all contributions to F{A(r)} with n gra- 
dients, one has to calculate terms of higher order than n in the iterative solution of the 
Gelfand-Dikii equation. In the case of complex A(r) KKSL derive a 2 x 2 matrix generaliza- 
tion of the Gelfand-Dikii equation, which is a form of the semiclassical Eilenberger equation 
0]. In this case no reshuffling of terms in the iterative procedure is necessary. To take care 
of a zero-mode KKSL only have to consider the next order in the gradient expansion. While 
in the work [Q gradients of the magnitude and the phase of A(r) are handled separately, an 
expansion in terms of gradients of A(r) itself is given in [^. 

In this work we shall develop an alternative algorithm for calculating the gradient ex- 
pansion of the free energy of a superconductor. We derive a pseudo-Schrodinger equation 
which is equivalent to the Eilenberger equation and develop a gradient expansion of the local 
density of states. In this expansion there appears a zero-mode which needs to be taken care 
of. In contrast to KKSL we fix this problem by making use of a non-linear constraint which 
naturally appears in our formalism. We do believe that this is an interesting and efficient 



alternative to the methods developed by KKSL. Our method works equally well for real and 
complex A(r) and directly generates the gradient expansion of the local density of states of 
the superconductor, from which the free energy can be obtained by simple integration. As 
an application of our method, we have explicitly calculated the local density of states p(r, u) 
and the free energy F{A(r)} up to fourth order in the gradients. We confirm the result by 
KKSL [^,^ for the fourth order correction to F{A(r)}, who obtained a discrepancy with 
the expression published by Tewordt B. 



II. FROM THE GORKOV EQUATIONS TO A PSEUDO-SCHRODINGER 

EQUATION 

In this section we shall map the problem of calculating the local density of states (DOS) 
of a superconductor onto the problem of calculating the time-evolution of the wave-function 
of an effective J = 1 quantum spin in a time- dependent complex magnetic field. For real A 
we have recently used a similar mapping to study the effect of order parameter fluctuations 



in disordered Peierls chains |10]. Here we generalize this mapping to the case of complex A. 



A. Reduction to an effective one-dimensional problem 

Starting point of our manipulations are the Gorkov equations for the matrix Green's 
function of a clean superconductor in the magnetic field H(r) = Vr x A(r), 

where o"o is the 2x2 unit matrix, and the differential operator Tir is given by 

[-iVr + -A(r)]2 

We have set ^ = 1 and measure energies and frequencies with respect to the chemical 
potential /i. Here c is the speed of light and — e is the charge of an electron. Note that KKSL 
0J^ start from the Bogoliubov - de Gennes equations (which are eigenvalue equations for 
the wave-functions of the superconductor) while our starting point are the Gorkov equations 
(which are differential equations for the Green's functions) . 

If |A(r)| ^ yU and if we are only interested in low-energy, long wavelength properties of 
the superconductor, we may linearize the energy dispersion. For calculating the local DOS 
and the free energy, we only need to know the Green's function at coinciding points r = r'. 
Following Waxman [|ll|, we write 



g{r, r, a;) ^ — (^„(r, r, u))^ . (3) 

Here u^ is the rf-dimensional DOS of free fermions at the Fermi energy (including both spin 
directions), and {■ ■ ■)^ = J ■ ■ ■ ^^ denotes directional averaging over the directions of the 
three-dimensional unit vector n. The auxiliary Green's function ^n(r, r','^) satisfies 



f uj- V{r) + ivpn ■ Vr A(r) \ 

l^ A*(r) u~V{r)-ivFn-Vr ) 

x^„(r,r',cu) = (5(n-(r-r')K, (4) 

where V{r) = -n ■ A(r) and vp is the Fermi velocity. Note that Eq.(D is the Green's 
function analog of the Andreev equations for the wave- function . Introducing r = xn + r^ 
and dx = n ■ Vr (where n ■ r^ = 0), dropping the dependence on the parameters r^ and n, 
and setting for simplicity vp = 1, we can write Eq.(^) as a truly one- dimensional equation, 

( uj-V{x)+id^ A(x) \ , , . 

[ A*{x) ij-V{x)-id,J^^'''^''^' 

= 6{x — x')cro • (5) 

To eliminate the forward scattering potential V{x) in Eq.(|D we set 

g{x, x', uj) = e-^^(")"='^i(x, x', c^)e^^("')'^3 , (6) 

where a^ is the usual Pauli matrix. Choosing the real function 6{x) such that 

^dxOix) = V{x) , (7) 



we obtain 



with 



A^xt -^-1 ) ^' ^'^' '^'' "^ = ^^"^ " "^'^^^ ' ^^^ 



A(x) = A(x)e'^(") . (9) 



B. Pseudo-Schrodinger equation for the local DOS 

The gradient expansion of the free energy can be obtained from the local DOS. Because 
Qi{x,x,u) has the same diagonal elements as Q{x,x,u), the local DOS pi(x,u;) of our 
effective one-dimensional model can be written as 

Pi{x,uj) = —'K^^lmTr[Qi{x,x,uj + iO)] , (10) 

where due to the trace we get a factor of 2 which takes both spin directions into account. Note 
that in the matrix equation (§) the derivative operator dx is proportional to the Pauli matrix 
(Ts. In the following it will be advantageous if the derivative operator dx is proportional to 
the unit matrix. Using cr| = (Tq, we obtain from Eq.(^ a differential equation with this 
property by setting [|lO| fl, 

g2{x,x',u) = a3Qi{x,x',uj) , (11) 



so that Eq.(|) implies 



= b{x — x')cro 



(12) 



For simplicity we have omitted the frequency label. The three Pauli matrices are denoted 
by Gi, 2 = 1, 2, 3, and a± = |[o"i ± ia2]- We now generalize the method described for real A 



in Ref. [|I0| to the case of complex A. We start by making the non-Ahelian Schwinger-ansatz 

SHI, 



Q2{x,x') = U{x)Qo{x,x')U (x) 



where U{x) is an invertible 2x2 matrix. It is easy to see that the solution of Eq.(12 
indeed be written in this form provided Qo and U satisfy 

[id^ao + ioas] Qo{x, x') = 6{x - x')cro , 



(13) 



can 



(14) 



idJJ{x) = uj[U{x)a3 - a3U{x)] 

+ [A{x)a+ - A*{x)a_]U{x) . 



(15) 



We parameterize U{x) as follows [|1^,0, 

This leads to the following system of equations for the complex functions $±(x) and $3(x), 



(9^$+ = -2icj$+ + A* - A<l>^ , 
d^<^^ = 2iuj<^_ - A[l - 2$+$_] 
(9^$3 = -iA$+ . 



(17a) 
(17b) 
(17c) 



For real A these equations reduce to the set of equations given in Ref. |]10[- Recently 
Schopohl [|1^ used an analogous procedure to map the semiclassical Eilenberger equations 
P] onto a similar set of non-linear equations. 

The one-dimensional local DOS defined in Eq.([TO|) can now be written as 

Pi{x, uj) = z/iRei?(x, uj + iO) , (18) 

where the complex variable R{x) is defined by (suppressing again the frequency label) 

i?(a;) = 1-2$+ (a;) <l>_(x) . (19) 



In principle one could now try to solve Eqs.( |17a] ) and ( |17bD iteratively in powers of the 
derivatives of A(a;), and then obtain the gradient expansion of the local DOS from Eqs.(|18|) 
and (0) . However, the structure of the iterative solution becomes more transparent if we 
introduce the three-component complex vector 



4'(x) 



Z+(x) 
R{x) 
Z(x) 



-v^[l-<l>+(a;)<l>„(a;)]<l>+(a;) 
l-2<l>+(s)<l>_(a;) 



(20) 



such that the 2x2 matrix Green's function at coinciding space points, 



is equal to 



g{x) = -[g2{x + o+,x) + g2{x,x + o+)] 



i 1 

g(x) = -R{x)a3 + -j=[Z_{x)a+ + Z+{x)a_ 



(21) 



(22) 



Because the components of ip can be parameterized by only two variables <l>±, they satisfy 
a constraint. Introducing the row vector 



i>'"ix) = [-Z4x),Rix),-Z4x)] , 
the constraint can be simply written as 

ip'^{x)ip{x) = R^{x) - 2Z+{x)Z_{x) = 1 
which can also be expressed as 

Differentiating the vector iIj{x) we find the equation of motion 

- d^ipi^x) = H{x)i){x) , 

with the pseudo-Hamiltonian given by 

/ 2iuj 72 A* (x) 

H{x) = V2A(x) V2A*(a;) 

V v^A(a;) -2iu 

Note that H{x) can also be written as 

H{x) = 2iujJ^ + A(x) J_ + A*(x) J+ 

where the Jj are spin J = 1 operators in the representation 



(23) 



(24) 



(25) 



(26) 



(27) 



(28) 



J:^ 



/I 


Vo -1 



J+ = V2 



/O 1 

1 1 = j! 
yo 



(29) 



Using the alternative parameterization oiip{x) in terms of the 2x2 matix g{x) given in (^ 
we find that our pseudo-Schrodinger equation is equivalent to the Eilenberger equation, 
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dxg{x) = [iucTs — 'iA(x)(T+ + iA*{x)a^,g{x)] . (30) 

This can be seen by using the commutation relations for the Pauh matrices and reducing 
Eq. (pOD to the equation of motion (pG]). Thus we have found an unconventional derivation 
of the Eilenberger equations which is based on the non-Abelian Schwinger-ansatz. 

Formally Eq.(P^) looks like the imaginary time Schrodinger equation for a J = 1 quantum 
spin subject to an imaginary time-dependent magnetic field. Recall that the real part of 
the second component of the state ip{x) can be identified with the local DOS. Because 
our pseudo- Schrodinger equation is linear, the gradient expansion of the local DOS can 
now be generated by a straightforward iterative calculation of the state iIj{x) in powers of 
gradients. Let us emphasize that, apart from the semiclassical approximation (Eq.(^), so 
far no approximation has been made. We have simply mapped the original problem onto 
an effective pseudo-Schrodinger equation, which is the most convenient starting point for 
setting up the gradient expansion. 

III. RECURSIVE ALGORITHM AND GRADIENT EXPANSION 

The gradient expansion of the local DOS is directly obtained from the second component 
of the gradient expansion oiiplx) . For convenience we develop the gradient expansion oiiplx) 
for imaginary frequencies u = lE, because then our pseudo-Hamiltonian ( pT]) is Hermitian 
and left and right eigenvectors are identical. Suppose we expand the solution of Eq. (^Hf) in 
the form 

oo 

Hx) = Y.M^), (31) 

n=0 

— * 

where by definition ipn{x) involves n derivatives with respect to x. Obviously 

H{x)i^o{x) = , (32) 

— * 

i.e. iI)q{x) must be an eigenvector of H{x) with eigenvalue zero. The existence of such an 
eigenvector follows trivially from the fact that our pseudo-Hamiltonian (^) can be inter- 
preted as the Zeeman-Hamiltonian of a J = 1 quantum spin in an external magnetic field. 
Note that Eq. (|32D determines iPq{x) only up to an overall multiplicative factor, which is 
fixed by requiring that the components of iIjq{x) satisfy the constraint (0). This yields 
(with oj = IE) 



^o(a;) 



^2 + |A( 



X] 



E 



(33) 



For the higher order terms we obtain the simple recursion relation 

d^ipnix) = -H{x)i)n+i{x) , n = 0,l,.... (34) 



Because one of the eigenvalues of H{x) vanishes, the inverse of H{x) does not exist, so that 
we cannot simply solve Eq.(|3^ by multiplying both sides by H~^{x). As a consequence, 
Eq.(^4D determines ipn+i{x) only up to a vector proportional to ^'0(2^)5 



1pn+l{x) = -H^ {x)d^i)n{x) + Cn+l{x)%{x) 



(35) 



where H^^{x) is the inverse of H{x) in the subspace orthogonal to il)o{x). Using the fact 
that the two non- vanishing eigenvalues of H{x) are given by ±2[£'^ + |A(a;)p]^/^, we find 



mh 



H(x) 



X] 



4[E2 + |A(a;) 



(36) 



i.e. Hj_^{x) is proportional to H{x). To fix the constant c„+i(a;) in Eq. (|35|) , we require that 
the components oi Y17=o i^ii^) satisfy the constraint (plf). This implies 



1 " ~ 

Cn+lix) = --J2'iplix)lpn+l-iix) 



(37) 



where the vector ipii^) is obtained from ipii^) by exchanging the first and third components 
and multiplying them by —1, see Eq.(^). For odd n we can show that Cn{x) = 0. We 
thus obtain an explicit and very compact recursive algorithm for calculating the gradient 
expansion of the local DOS. To zeroth order the vector ?/'(x) is given by Eq. (P^. This 
corresponds to the adiabatic approximation of elementary quantum mechanics. The step 
ra — > ra + 1 is summarized as follows 



ipi given for i = 0, . . . ,n 



(38) 



It is easy to implement this iterative algorithm on a symbolic manipulation program (such 
as Mathematica) . In this way the lowest few terms in the gradient expansion can be obtained 
in a straightforward manner. 

Given the gradient expansion of the local DOS (which can be directly obtained from the 
second component of ip{x)), we can calculate the free energy by simple integrations. Using 
the fact that in the normal state R = 1, the difference between the free energy densities 
in the superconducting and normal state of our effective one-dimensional model at inverse 
temperature (3 is given by uiflx), where 



/(^) 



-Re 



Re 



/ dtu [R{x, to + iO) -I] \n(l + e 
00 ]^ 



-I3u 



271 






(39) 



Here Un = (2n + l)7r//5 are fermionic Matsubara frequencies, and 



N{x,uj) 



duj'R{x, uo') . 



(40) 



Note that z/iRe[A^(a;, a; + ^0)] is the integrated one- dimensional local DOS. In Eq.(p9D we 
have used the fact that i?(x, u) is analytic in the upper half of the complex cu-plane. The 
local DOS of the three-dimensional superconductor is then given by 



p(r, u) = z/3 (Rei?n(r, uj + iO))^ 



(41) 



where we have now written -Rn(r, uj) for R{x, u) to indicate the dependence on all parameters. 
Similarly we write /n(r) for f{x) and obtain for the difference between the free energies of 
the three dimensional system in the superconducting and normal state, 



F{A(r)}= fdh 



1^3 (/n(r))„ + 



lAfr) 



A 



+ 



[H(r) - He(r)]^ 



ill 



(42) 



The second term is the field energy of the superconducting pair potential (where A is the 
coupling constant of the Gorkov electron-electron interaction), and the last term is the 
magnetic field energy of the superconductor (where He(r) is the externally applied magnetic 
field i). 

The above equations allow for a simple recursive calculation of the gradient expansion of 
the free energy. To compare our results with Tewordt and KKSL 0,^ we use Eq.(^) to 
calculate all terms in the gradient expansion up to fourth order. Systematically adding total 
derivatives to the expressions for the free energy (which do not change the bulk properties 
of the superconductor), we find the following expressions: 



Zeroth order: 



p^'^{r,uj) = use{uj' 



lAI 



u; 



c<;2-|A|2 



(43) 



27r ^ 

f^ UJ„>0 



F(o){A(r)} = / dh 

Second order: 

p^^\r,uj) = Us\uj\e{uj^-\A\^ 



cu2 + |A|2-cJ„ 



+ 



|A(r)p , [H(r) - H,(r)]^ 



A 



in 



■ (44) 



5 [n ■ Vr|A|2]^ 1 [(n ■ Vr)^\A\^] - 3|n ■ DrA| 



32fcj2- lAP' 



fcj2-|A|2)2 



(45) 



F(^){A(r)} = .3/rf^r^E 



i^n>0 



1 [n-VJA 



21 2 



32 



1 |n-D.A| 



;.;2 + |A|2)^ 8(^2 + |A|2). 



(46) 



Fourth order: 

\2048(^2„|A|2)- 512 (cc;2-|A|2)- 

7 5|n • DrA|4 + 4 [n • Vr|A|2] [(n • Vr)3|A|2] - 10 [n • Vr (|n • DrAp [n • Vr|A|2])] + 3 [(n • Vr)^\A\^y 



128 
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(w2_|A|2)2 



^ 1 [(n ■ V.)^|A|2] - 5 [(n ■ V,)2|n ■ D.Ap] + 5|(n ■ D,)2A|2 \ 
32 (c.2-|A|2)5 // 

.(^) {A(r)} . .3 / .^rf E (-^^^^^^ - ^ K^,M!^V^ 
J /3^^n\ 2048^^2 ,|A|2^- 512 f^,? , |a|2)2 



w„>0 



K + |A|2)^ Oi^ K + |A|2)3 



^ 1 5|n.D,A|4-10|n.D,A|2[(n.V,)2|A|2] + [(n.V,)2|A|2]^ 1 |(n ■ D,)2A|2 \ 
128 K + |A|2)i 32(^2 + |A|2)l// 

Here Q{x) is the step function, and Dp = Vr + 2z-A(r) is the covariant gradient. In deriving 
the above expressions we have used 

|(n-Vr)'A| = |(n-Dr)'A| , / = 0,1,2,... . (49) 

Note that all terms in our expansion are gauge invariant. Since the terms involving an odd 
number of derivatives cancel after directional averaging we only presented the terms of even 
order. The directional averaging is easily done using 

((n-A)(n-B))„ = lA-B, (50) 

((n . A)(n • B)(n ■ C)(n • D))„ = ^ [(A • B) (C ■ D) 

+ (A • C) (B ■ D) + (A • D) (B ■ C)] . (51) 

Because after directional averaging the above expressions look even more complicated, we 
do not give the explicit results in this work. Comparing Eq.(^) with the fourth order 
correction to the free energy given by Tewordt and by KKSL @J^, we find that our 
result agrees with that of KKSL (up to a total derivative which eliminates the first term on 
the right-hand-side of Eq. (0)). We have not been able to find in the published literature 
an explicit expression for the fourth order correction to the gauge invariant local DOS, so 
that we cannot compare Eq.(^) with the results of other authors. 

IV. CONCLUSION 

We have developed an efficient algorithm for calculating the gradient expansion of the 
local DOS and the free energy of a superconductor in an external magnetic field. The 
linearization of the energy dispersion allows to obtain the local DOS directly from the 
second component of a three component wave-function ip^x) satisfying a pseudo-Schrodinger 
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equation, which is formally identical to the Schrodinger equation of a J = 1 quantum 
spin in a time- dependent complex magnetic field. The gradient-expansion of this pseudo- 
Schrodinger equation turns out to be quite simple and straightforward and can be easily 
implemented on a symbolic manipulation program such as Mathematica. It seems to us 
that our algorithm gives an interesting alternative to the algorithms developed by KKSL 
0J^. Our final result for the contribution to the gauge-invariant free energy containing four 
gradients agrees with KKSL [Q and is in disagreement with the expression given by Tewordt 
0. We believe that our algorithm could also prove advantageous in other problems where 
fluctuation effects can be mapped onto a 2 x 2 matrix equation for a Green's function in 
external fields. One interesting example are two coupled Tomonaga-Luttinger models. 

We would like to thank Simon Kos for pointing out to us that our pseudo-Schrodinger 
equation is in fact equivalent to the Eilenberger equation. This work was financially sup- 
ported by the DFG (Grants No. Ko 1442/3-1 and Ko 1442/4-1). 
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